Exploration

Packages needed

library(naniar) # for missing data viz
library(tidyverse) # for data wrangling
library(skimr) # for statistical summaries
library(DT) # for nice tables 
library(mice) # for missing data imputation

set.seed(123) # for reproducibility

Riskfactors dataset

The data is a subset of the 2009 survey from BRFSS (Behavioral Risk Factor Surveillance System), an ongoing data collection program designed to measure behavioral risk factors for the adult population (18 years of age or older) living in households.

More information: https://www.rdocumentation.org/packages/naniar/versions/0.5.1/topics/riskfactors

skim(riskfactors)
Data summary
Name riskfactors
Number of rows 245
Number of columns 34
_______________________
Column type frequency:
factor 18
numeric 16
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
state 0 1.00 FALSE 52 6: 12, 34: 12, 25: 11, 12: 10
sex 0 1.00 FALSE 2 Fem: 153, Mal: 92
marital 1 1.00 FALSE 6 Mar: 131, Wid: 41, Div: 39, Nev: 24
pregnant 215 0.12 FALSE 2 No: 28, Yes: 2
education 1 1.00 FALSE 6 4: 85, 6: 75, 5: 59, 3: 14
employment 0 1.00 FALSE 7 1: 92, 7: 80, 2: 21, 8: 18
income 0 1.00 FALSE 10 >75: 46, 35-: 38, 50-: 32, 25-: 28
veteran 3 0.99 FALSE 5 5: 198, 3: 35, 4: 6, 2: 2
hispanic 2 0.99 FALSE 2 No: 221, Yes: 22
health_general 0 1.00 FALSE 6 Goo: 82, Ver: 56, Exc: 45, Fai: 40
health_cover 0 1.00 FALSE 2 Yes: 217, No: 28
provide_care 3 0.99 FALSE 2 No: 181, Yes: 61
activity_limited 3 0.99 FALSE 2 No: 170, Yes: 72
drink_any 2 0.99 FALSE 2 No: 132, Yes: 111
smoke_100 2 0.99 FALSE 2 No: 126, Yes: 117
smoke_days 128 0.48 FALSE 3 Not: 84, Eve: 26, Som: 7
smoke_stop 212 0.13 FALSE 2 Yes: 19, No: 14
smoke_last 161 0.34 FALSE 6 7: 59, 5: 9, 6: 6, 8: 5

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 0 1.00 58.11 17.50 7.00 48.00 59.00 70.00 97.00 ▁▃▇▇▃
weight_lbs 10 0.96 174.27 44.70 96.00 144.00 170.00 195.00 410.00 ▆▇▂▁▁
height_inch 2 0.99 66.35 3.80 57.00 64.00 66.00 69.00 75.00 ▂▇▇▇▃
bmi 11 0.96 27.78 6.57 14.21 23.35 26.66 31.32 55.72 ▃▇▃▁▁
children 0 1.00 0.42 0.88 0.00 0.00 0.00 0.00 4.00 ▇▁▁▁▁
health_physical 0 1.00 4.12 9.11 -9.00 0.00 0.00 3.00 30.00 ▁▇▁▁▁
health_mental 0 1.00 3.14 8.19 -9.00 0.00 0.00 1.00 30.00 ▁▇▁▁▁
health_poor 113 0.54 5.46 10.76 -9.00 0.00 0.00 5.00 30.00 ▁▇▁▁▂
drink_days 134 0.45 9.29 10.20 -9.00 2.00 4.00 15.00 30.00 ▁▇▃▃▂
drink_average 135 0.45 1.41 2.57 -9.00 1.00 1.00 2.00 9.00 ▁▁▇▆▁
diet_fruit 8 0.97 324.71 338.11 -9.00 60.00 209.00 365.00 2555.00 ▇▁▁▁▁
diet_salad 8 0.97 175.48 140.00 -9.00 52.00 156.00 261.00 730.00 ▇▇▃▁▁
diet_potato 8 0.97 111.68 132.42 -9.00 24.00 104.00 156.00 1095.00 ▇▁▁▁▁
diet_carrot 8 0.97 91.52 101.32 -9.00 24.00 52.00 156.00 365.00 ▇▂▂▁▁
diet_vegetable 8 0.97 407.15 316.25 -9.00 156.00 365.00 626.00 1825.00 ▇▇▃▁▁
diet_juice 8 0.97 227.85 307.28 -9.00 24.00 104.00 365.00 2555.00 ▇▁▁▁▁ 
n_var_miss(riskfactors)
## [1] 24

Missingess patterns

gg_miss_upset(riskfactors)

The documentation on {naniar} says that the default option of gg_miss_upset is taken from UpSetR::upset - which is to use up to 5 sets and up to 40 interactions.

Summaries

Visual summaries

gg_miss_var(riskfactors)

gg_miss_var(riskfactors, show_pct = TRUE)

gg_miss_var(riskfactors,
            show_pct = T,
            facet = sex)

gg_miss_fct(x = riskfactors, fct = marital)

Data summaries

riskfactors %>%
  group_by(marital) %>%
  miss_var_summary() %>% 
  DT::datatable()

MICE

Motivation

We want to assess the predictors of health_physical, a continuous variable corresponding to the number of days in the last 30 days for which the respondent said their physical health was good.

Explanatory variables of interest include the ones related to diet as well as health_poor. If we were to run a Complete Case Analysis regression analysis, we’d effectively be throwing away nearly half our sample.

Preprocessing

riskfactors <- riskfactors %>% select(-c("pregnant", starts_with("smoke"), "drink_days", "drink_average"))

Impute

  • To perform imputations, you pass the dataset name to mice(). The default number of imputations is m = 5.
  • the quickpred() can be passed on to the pred argument if the researcher is interested in quickly selecting predictors in datasets that contain many variables. The threshold for correlation can be specified with mincor.

Here we select predictors with a minimum correlation of \(\rho=.30\) :

impute <- mice(riskfactors, 
               m = 5, 
               pred = quickpred(riskfactors, mincor = .3),
               print = FALSE)

Checking the class of impute

class(impute)
## [1] "mids"

Check the convergence

plot(impute)

The plot shows the mean (left) and standard deviation (right) of the imputed values only. In general, we would like the streams to intermingle and be free of any trends at the later iterations.

Check the method used

As the documentation shows (https://www.rdocumentation.org/packages/mice/versions/3.5.0/topics/mice), many methods can be specified for imputation. To check which one our algorithm chose:

cbind(impute$meth)
##                  [,1]     
## state            ""       
## sex              ""       
## age              ""       
## weight_lbs       "pmm"    
## height_inch      "pmm"    
## bmi              "pmm"    
## marital          "polyreg"
## children         ""       
## education        "polyreg"
## employment       ""       
## income           ""       
## veteran          "polyreg"
## hispanic         "logreg" 
## health_general   ""       
## health_physical  ""       
## health_mental    ""       
## health_poor      "pmm"    
## health_cover     ""       
## provide_care     "logreg" 
## activity_limited "logreg" 
## drink_any        "logreg" 
## diet_fruit       "pmm"    
## diet_salad       "pmm"    
## diet_potato      "pmm"    
## diet_carrot      "pmm"    
## diet_vegetable   "pmm"    
## diet_juice       "pmm"

Looks like the algorithm used predictive mean matching (pmm) for the numerical variables, and either polytomous regression (polyreg) for unordered categorical variables or logistic regression for the binary categorical variables (logreg)

We can also change the method. Let’s use Bayesian linear regression - norm - for the bmi variable:

meth <- impute$meth
meth["bmi"] <- "norm"
cbind(meth)
##                  meth     
## state            ""       
## sex              ""       
## age              ""       
## weight_lbs       "pmm"    
## height_inch      "pmm"    
## bmi              "norm"   
## marital          "polyreg"
## children         ""       
## education        "polyreg"
## employment       ""       
## income           ""       
## veteran          "polyreg"
## hispanic         "logreg" 
## health_general   ""       
## health_physical  ""       
## health_mental    ""       
## health_poor      "pmm"    
## health_cover     ""       
## provide_care     "logreg" 
## activity_limited "logreg" 
## drink_any        "logreg" 
## diet_fruit       "pmm"    
## diet_salad       "pmm"    
## diet_potato      "pmm"    
## diet_carrot      "pmm"    
## diet_vegetable   "pmm"    
## diet_juice       "pmm"

We must re-run the imputation:

impute <- mice(riskfactors, meth = meth, print = FALSE)

We can plot the trace lines again to see what the convergence looks like:

plot(impute)

Increase the number of imputations

We’ve seen above that the lines do intermingle well, especially at the later iterations with just 5 iterations. Sometimes, it is useful to increase the number of iterations just to confirm that there is indeed no trend. We can increase the number of iterations to 20 by running 15 additional iterations using the mice.mids() function:

impute_20 <- mice.mids(impute, maxit = 15, print = FALSE)
plot(impute_20)

Further diagnostic testing

We can check our imputations for certain variables by comparing them to observed values - under the MCAR assumption, the imputations should indeed have the same distribution as the observed data. Under MAR, the distributions may be different but nevertheless very large differences warrant further investigation.

  • For health_poor days:
stripplot(impute, health_poor ~ .imp, pch = 20, cex = 2)

  • For bmi:
stripplot(impute, bmi ~ .imp, pch = 20, cex = 2)

bmi was imputed with Bayesian linear regression and (the range of) imputed values looks a little different than the observed values but there is still overlap.

Perform regression analysis and pool estimates

Finally, we need to run the regression on each of the 5 datasets and pool the estimates together to get average regression coefficients and correct standard errors. The with() function in the mice package allows us to do this.

fit <- with(impute, lm(health_physical ~ age + bmi + sex + marital + children + employment + health_mental + activity_limited + provide_care + health_poor + health_general + diet_vegetable + diet_salad + diet_potato + diet_juice + diet_fruit + diet_carrot))

fit
## call :
## with.mids(data = impute, expr = lm(health_physical ~ age + bmi + 
##     sex + marital + children + employment + health_mental + activity_limited + 
##     provide_care + health_poor + health_general + diet_vegetable + 
##     diet_salad + diet_potato + diet_juice + diet_fruit + diet_carrot))
## 
## call1 :
## mice(data = riskfactors, method = meth, printFlag = FALSE)
## 
## nmis :
##            state              sex              age       weight_lbs 
##                0                0                0               10 
##      height_inch              bmi          marital         children 
##                2               11                1                0 
##        education       employment           income          veteran 
##                1                0                0                3 
##         hispanic   health_general  health_physical    health_mental 
##                2                0                0                0 
##      health_poor     health_cover     provide_care activity_limited 
##              113                0                3                3 
##        drink_any       diet_fruit       diet_salad      diet_potato 
##                2                8                8                8 
##      diet_carrot   diet_vegetable       diet_juice 
##                8                8                8 
## 
## analyses :
## [[1]]
## 
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children + 
##     employment + health_mental + activity_limited + provide_care + 
##     health_poor + health_general + diet_vegetable + diet_salad + 
##     diet_potato + diet_juice + diet_fruit + diet_carrot)
## 
## Coefficients:
##            (Intercept)                     age                     bmi  
##              1.739e+00              -6.337e-02               1.503e-01  
##              sexFemale         maritalDivorced          maritalWidowed  
##              2.048e+00               1.586e-01               1.593e+00  
##       maritalSeparated     maritalNeverMarried  maritalUnmarriedCouple  
##             -4.300e+00              -5.817e-01              -3.503e+00  
##               children             employment2             employment3  
##             -6.484e-01               5.212e-01              -2.866e+00  
##            employment4             employment5             employment7  
##             -2.320e+00               1.163e+00              -4.259e-01  
##            employment8           health_mental      activity_limitedNo  
##             -2.402e+00               1.047e-01              -4.212e+00  
##         provide_careNo             health_poor  health_generalVeryGood  
##              4.186e-01               1.132e-01               4.179e-01  
##     health_generalGood      health_generalFair      health_generalPoor  
##              6.289e-01               3.133e+00               1.679e+01  
##  health_generalRefused          diet_vegetable              diet_salad  
##             -8.066e+00               5.988e-04               1.030e-03  
##            diet_potato              diet_juice              diet_fruit  
##             -5.565e-03               4.084e-03              -4.019e-05  
##            diet_carrot  
##             -1.476e-03  
## 
## 
## [[2]]
## 
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children + 
##     employment + health_mental + activity_limited + provide_care + 
##     health_poor + health_general + diet_vegetable + diet_salad + 
##     diet_potato + diet_juice + diet_fruit + diet_carrot)
## 
## Coefficients:
##            (Intercept)                     age                     bmi  
##              0.9990222              -0.0686284               0.1778563  
##              sexFemale         maritalDivorced          maritalWidowed  
##              1.7289946               0.1417556               1.5580613  
##       maritalSeparated     maritalNeverMarried  maritalUnmarriedCouple  
##             -2.7448279              -0.4858848              -3.5094290  
##               children             employment2             employment3  
##             -0.6290995               1.0762904              -2.5990223  
##            employment4             employment5             employment7  
##             -1.6601152               0.9011993              -0.8981470  
##            employment8           health_mental      activity_limitedNo  
##             -2.9122112               0.1028946              -4.4034850  
##         provide_careNo             health_poor  health_generalVeryGood  
##              0.2056059               0.1056196               1.1528270  
##     health_generalGood      health_generalFair      health_generalPoor  
##              1.1309828               3.6589914              17.4273975  
##  health_generalRefused          diet_vegetable              diet_salad  
##             -6.9687402              -0.0000866               0.0043914  
##            diet_potato              diet_juice              diet_fruit  
##             -0.0041553               0.0038780              -0.0001935  
##            diet_carrot  
##             -0.0011747  
## 
## 
## [[3]]
## 
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children + 
##     employment + health_mental + activity_limited + provide_care + 
##     health_poor + health_general + diet_vegetable + diet_salad + 
##     diet_potato + diet_juice + diet_fruit + diet_carrot)
## 
## Coefficients:
##            (Intercept)                     age                     bmi  
##               2.704129               -0.071143                0.135268  
##              sexFemale         maritalDivorced          maritalWidowed  
##               1.830159                0.200566                1.503719  
##       maritalSeparated     maritalNeverMarried  maritalUnmarriedCouple  
##              -4.168948               -0.415450               -3.403838  
##               children             employment2             employment3  
##              -0.749504                0.971755               -2.642386  
##            employment4             employment5             employment7  
##              -2.136584                0.841615               -0.552605  
##            employment8           health_mental      activity_limitedNo  
##              -2.480937                0.104532               -4.453114  
##         provide_careNo             health_poor  health_generalVeryGood  
##              -0.030878                0.106614                0.647872  
##     health_generalGood      health_generalFair      health_generalPoor  
##               0.988748                3.445597               16.724057  
##  health_generalRefused          diet_vegetable              diet_salad  
##              -7.735275                0.001092                0.002808  
##            diet_potato              diet_juice              diet_fruit  
##              -0.003753                0.004438               -0.001017  
##            diet_carrot  
##              -0.002032  
## 
## 
## [[4]]
## 
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children + 
##     employment + health_mental + activity_limited + provide_care + 
##     health_poor + health_general + diet_vegetable + diet_salad + 
##     diet_potato + diet_juice + diet_fruit + diet_carrot)
## 
## Coefficients:
##            (Intercept)                     age                     bmi  
##              1.5677271              -0.0764285               0.1665056  
##              sexFemale         maritalDivorced          maritalWidowed  
##              1.4635227               0.3138129               1.8212133  
##       maritalSeparated     maritalNeverMarried  maritalUnmarriedCouple  
##             -3.0856733              -0.7627401              -3.1658295  
##               children             employment2             employment3  
##             -0.7966279               0.1802052              -2.6618967  
##            employment4             employment5             employment7  
##             -1.2421182               0.2962280              -1.0384293  
##            employment8           health_mental      activity_limitedNo  
##             -2.8780833               0.0864878              -3.8797074  
##         provide_careNo             health_poor  health_generalVeryGood  
##              0.3103232               0.1812182               1.3255393  
##     health_generalGood      health_generalFair      health_generalPoor  
##              1.1020951               4.0242406              16.8904464  
##  health_generalRefused          diet_vegetable              diet_salad  
##             -5.9716591               0.0006650               0.0050499  
##            diet_potato              diet_juice              diet_fruit  
##             -0.0071418               0.0032890              -0.0006714  
##            diet_carrot  
##              0.0009036  
## 
## 
## [[5]]
## 
## Call:
## lm(formula = health_physical ~ age + bmi + sex + marital + children + 
##     employment + health_mental + activity_limited + provide_care + 
##     health_poor + health_general + diet_vegetable + diet_salad + 
##     diet_potato + diet_juice + diet_fruit + diet_carrot)
## 
## Coefficients:
##            (Intercept)                     age                     bmi  
##              2.6453979              -0.0764877               0.0868773  
##              sexFemale         maritalDivorced          maritalWidowed  
##              1.8208563               0.3403289               2.0334635  
##       maritalSeparated     maritalNeverMarried  maritalUnmarriedCouple  
##             -6.8832296              -1.2065212              -4.8079140  
##               children             employment2             employment3  
##             -0.6356130               0.7207182              -2.9229650  
##            employment4             employment5             employment7  
##             -2.0566007               0.3248994              -0.9867242  
##            employment8           health_mental      activity_limitedNo  
##             -2.6273493               0.0906184              -3.4591495  
##         provide_careNo             health_poor  health_generalVeryGood  
##              0.4729385               0.2205807               1.4180455  
##     health_generalGood      health_generalFair      health_generalPoor  
##              2.0186421               4.4735674              16.6651833  
##  health_generalRefused          diet_vegetable              diet_salad  
##             -5.3054085               0.0007200               0.0037058  
##            diet_potato              diet_juice              diet_fruit  
##             -0.0071866               0.0039170              -0.0005845  
##            diet_carrot  
##              0.0004970

The fit object contains the regression summaries for each data set. The new object fit is actually of class mira (multiply imputed repeated analyses). We can double-check:

class(fit)
## [1] "mira"   "matrix"

Pool the analyses from object fit

This is the last step, where we pool the estimates from all 5 complete datasets to get average regression coefficients and correct standard errors.

pool.fit <- pool(fit)

summary_poolfit <- as.data.frame(summary(pool.fit))

summary_poolfit %>% DT::datatable()
class(pool.fit)
## [1] "mipo"       "data.frame"